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SUMMARY 

An  equation  of  missile  motion  in  the  pitching  plane  was  found, 
assuming  a  non-linearity  in  aerodynamic  stiffness  in  the  form  of  a 
cubic  term.   An  iterative  method  for  approximating  the  angle-of- 
attack  response  to  a  step  displacement  of  elevator  angle  was 
formulated.   The  results  of  applications  of  this  method  were 
compared  with  computer  solutions;  convergence  of  the  approximations 
Me^e  demonstrated. 

A  simple  feedback  control  system  was  designed,  and  the  iterative 
method  was  extended  to  the  case  of  a  control  system  step  input. 

The  method  was  thought  to  have  limited  practical  value  because 
of  the  labor  required  in  using  it. 
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AN  ITERATIVE  METHOD  FOR  FINDING  THE 
TRANSIENT  RESPONSE  OF  M  AUTOMATICALLY 
CONTROLLED  MISSILE  HAVING  NON-LINEAR  STIFFNESS 

INTRODUCTION 


Important  non-linearities  which  are  often  present  in  low  aspect 
ratio  configurations  are  those  of  the  lift  versus  incidence  and  the 
pitching  moment  versus  incidence  characteristics.   In  missile  design 
these  non-linearities,  as  well  as  any  others  present,  are  likely  to 
be  ignored  initially,  simply  because  linear  system  design  is  well 
developed  and  relatively  easy.   Automatic  computing  techniques  may 
then  be  used  to  determine  the  behaviour  of  the  system  when  the  non- 
linearities  are  included. 

It  would  be  a  great  help  for  the  designer  to  be  able  to  find 
the  transient  response  of  a  contemplated  design  without  going  to  the 
computer.   It  was  for  this  reason  that  the  present  investigation 
was  undertaken.   By  expressing  the  pitching  moment  as  the  sum  of  a 
linear  term  and  a  cubic  term  in  angle-of -attack,  ignoring  the 
non-linearity  in  lift  because  it  is  likely  to  be  less  important, 
an  equation  amenable  to  iterative  solution  can  be  formulated.   The 
first  approximation  to  the  transient  response  may  be  expressed  as  a 
Fourier  series,  which  must  be  cubed  to  perform  the  first  iteration. 
The  main  effort  of  the  development  will  be  to  formulate  a  practical 
method  for  performing  the  cubing  process,  for  the  first  and 
subsequent  iterations,  and  to  systematise  the  steps  necessary  for 
carrying  out  the  iterations. 

First,  an  equation  relating  angle-of -at tack  to  elevator  angle 
will  be  derived,  and  the  response  to  a  step  displacement  of  elevator 


. 


angle  will  be  found.   Later,  a  simple  automatic  feedback  control 
system  will  be  designed  for  the  purpose  of  formulating  a  method 
of  finding  the  response  to  a  control  system  input. 

The  author  takes  this  opportunity  to  express  his  gratitude  for 
the  generous  assistance  given  by  Mr.  P.  A.  T.  Christopher,  who 
instigated  and  supervised  this  investigation. 


ENGLISH-LETTER  NOMENCLATURE 

a  coefficient  of  restoring  term  in  Sec.  3 

a,  coefficient  of  damping  term  in  Sec.  3 

a.  .,  coefficients  of  terms  in  the  cube  of  Fourier  series 

in  Sec.  2 

A  gain  of  frequency  response  at  frequency  n  ^  p 

b  coefficient  of  cubic  term  in  Sec.  3 

b  coefficients  in  the  Fourier  series  for  "^  1 0 

B  phase  shift  of  frequency  response  at  frequency  mdF 

c    coefficients  in  Sec.  2 
n 


C    phase  angles  in  the  Fourier  series  for   °<  i 


n 


v 


C(s)  1/Kp  times  the  transfer  function  of  the  compensator- 
amplifier 

CT  lift  coefficient 

CL.  pitching  moment  coefficient 

C-,  partial  derivative  of  II.  with  respect  to  <x? \j 

Cp  coefficient  of  vane  damping  term 

db  decibels 

D  d/dt 

D  missile  body  diameter 

E  actuating  error 

E  coefficients  in  trigonometric  identities 

F  phase  angles  in  trigonometric  identities 


H(s)  transfer  function  used  in  closed-loop  iteration 

i  f(-l) 

I  moment  of  inertia  of  missile  about  y  axis 

L,  moment  of  inertia  of  vane  about  its  e.g. 

K,  zero  frequency  gain  of  the  aerodynamic  transfer 

function 

Kp  zero  frequency  gain  of  compensator-amplifier 

IC,  gain  of  input  amplifier 

ly.  distance  between  missile  e.g.  and  vane  e.g. 

L  lift 

m  mass  of  missile 

M  aerodynamic  pitching  moment  on  missile  about  y  axis 

M„  aerodynamic  pitching  moment  on  vane  about  its  e.g. 

M'(w)  component  of  moment  due  to  w 

M"(o<)  component  of  moment  due  to  cK 

M^  coefficient  of  the  linear  term  in  the  cubic  expression 
for  pitching  moment 

M^?  coefficient  of  the  cubic  term  in  the  cubic  expression 
for  pitching  moment 

M-  partial  derivative  of  M  with  respect  to  * 

n  harmonic 

P(s)  l/K-j  times  aerodynamic  transfer  function 

q  angular  velocity  of  missile  about  the  y  axis 

Q  step  function  in  Sec.  3 


Q  stagnation  pressure  in  appendices  I  and  VI 

IQ(  magnitude  of  step  function 

Q  non-dimensional  step  input 

|Q  j  magnitude  of  non-dimensional  step  input 

E.  ratio  of  non-linear  to  linear  steady-state  values 

s  Laplace  transform  variable 

S  wing  area 

t  time 

t,  time  measured  from  start  of  second  half  cycle 

T  period  of  square  wave 

T(d)  differential  operator  in  equation  for  missile 
aerodynamics 

u  component  of  missile  e.g.  velocity  in  x  direction 

U(D)  differential  operator  in  equation  for  missile 
aerodynamics 

V  velocity  of  missile 

V(s)  vane  transfer  function 

w  component  of  missile  e.g.  velocity  in  z  direction 

W(s)  actuator  transfer  function 

Z  aerodynamic  force  component  in  the  z  direction 

Z  partial  derivative  of  Z  with  respect  to  w 

Z^  partial  derivative  of  Z  with  respect  to   <K 

Zft  partial  derivative  of  Z  with  respect  to   n 


GREEK-LETTER  NOMENCLATURE 

o<    missile  angle-of -attack 
c<^         demanded  missile  angle-of -at tack 


\a<  |   magnitude  of  input  step  function 


o^^  reference  input 

c^-  vane  angle-of -at  tack 

o<    o(      ....  approximations  to   °4 

©<*.©<.*  corrections  to  approximations 

o^,  response  of  linear  system  to  unit  step  input 

/3  coefficient  of  the  cubic  term 

y  angle  between  missile  axis  and  vane  axis 

3  damping  ratio 

t  damping  ratio  of  missile 
i 

£  damping  ratio  of  vane 

$  damping  ratio  of  elevator  actuator 

h  elevator  deflection  angle 

1 1 1  magnitude  of  step  elevator  input 

In  J  demanded  elevator  angle 

0  angle  between  missile  axis  and  space  reference 

Q  angle  between  vane  axis  and  space  reference 

Z  time  constant  in  the  missile  aerodynamics  transfer 
function 

t*>p  fundamental  frequency  of  square  wave 

lO  natural  frequency  of  linear  system 


U^    natural  frequency  of  linear  missile 
u)  natural  frequency  of  vane 


L 


u^,    natural  frequency  of  elevator  actuator 


SUBSCRIPTS  NOT  SHOW  ELSEWHERE 

a  evaluated  function 

b  function  to  be  derived  from  an  evaluated  function 

L  linear 

N  non-linear 

S  steady-state 
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SECTION  1:  THE  EQUATION  OF  MISSILE  MOTION  IN  PITCH 

The  system  of  axes  and  sign  conventions  to  be  used  are 
illustrated  in  Fig.  1.   Body  axes  are  used.   The  missile  to  be 
considered  has  a  fixed  wing  and  an  all-moving  tail.   It  may  or  may 
not  be  of  cruciform  configuration,  since  pure  pitching  motion  will 
be  assumed. 

The  rigid  body  equations  of  motion  in  the  z  direction  and  about 
the  y  axis  (which  passes  through  the  center  of  mass)  are, 
respectively, 

Z  =  m(w  -  qu)  (l-l) 

and 

M  =  Iq  (1-2) 

The  same  equations  would,  of  course,  arise  if  the  general  three 
dimensional  equations  of  motion  were  written  down  and  the  confining 
assumptions  of  the  present  motion  imposed. 

The  component  in  the  z  direction  of  the  total  force  on  the 
missile  is  taken  as 

Z  =  V  +  zi  n  U-3) 

The  simplifying  assumption  has  been  made  here  that  the  gravitational 
force  may  be  neglected  because  it  is  small  compared  with  the 
aerodynamic  forces  available.   The  contribution  to  Z  due  to  q  is 
omitted  on  the  assumption  that  it  is  small.   Its  inclusion  would 
not  change  the  nature  of  the  analysis,  but  it  would  be  dropped 
eventually  anyway  because  values  for  it,  for  the  configuration  and 
flight  conditions  considered,  are  not  available.   Also  implicit  in 
Eq.  1-3  is  the  assumption  that  the  malalignment  force  is  zero, 
i.e.  Z  is  zero  when  w  and  n    are  zero. 


The  aerodynamic  moment  is  taken  as 

M  =  M'(w)  +  M  n  (1-4) 

where  M'(w)  is  to  be  read,  "The  component  of  moment  due  to  w." 
This  is  the  non-linear  term.   The  same  remarks  as  those  following 
Eq.  1-3  apply  to  the  exclusion  of  components  of  moment  due  to  q 
and  to  malalignment. 

Now  we  take 


■  u  =  V 

which  is  valid  for  small  angles  of  attack,  and  we  assume  that 
V  is  constant.   Then,  for  small  angles  of  attack, 


and 


Z   = 

w 


11 
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Substituting  Eqs.  1-5  and  1-6  in  Eq.  1-1  gives 

Z  =  m(V<*  -  qV) 
Substituting  Eqs.  1-6  and  1-7  in  Eq.  1-3  gives 

Also  we  replace  M1  (w)  by  M"(«<),  where  M"(<x)  is  to  be  read, 
"The  component  of  moment  due  to  o<  ."   Thus  Eq.  1-4  becomes 

M  =  M"  (°<)  +  M  h 
Combining  Eqs.  1-8  and  1-9  gives 

Z^  <*  +  Z  n  =  m(V<*  -  qV) 
Combining  Eqs.  1-2  and  1-10  gives 

M»(o<)  +  M  ^  =  Iq 


(1-5) 


(1-6) 
(1-7) 

(1-8) 
(1-9) 


(1-10) 


(1-11) 


(1-12) 
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Eliminating  q  between  Eqs.  1-11  and  1-12  gives 

mlV  '<<   -  IZ^<*  -  MVM"(<*)  =  IZ  A  +  mVM^  »j         (l-lj) 

Taking 

M"(«*)  =  M^  <*  +  M^->  e<*  (1-14) 

equation  1-13  becomes 

mlV  ^<  -  IZ^  <*  -  mVM^  K  -  mVM^*  °<^  =  1ZX\     +  mVMi,  1    U-15) 

Dividing  Eq.  1-15  through  by  mlV  gives 

Z_    M     M~^      Z„    M 
<*   "  <X~  <*  -  **  =  h  +  __  n  (1-16) 

mV     I      I       mV  *    I 
The  constants   ^  ,    $,  ,  K,,   £"  and  /3  are  now  defined 

(1-17) 
(1-18) 
(1-19) 
(1-20) 

(1-21) 

I 

Making  these  substitutions  in  Eq.  1-16  gives 

o<  +.  z  r,  u/^(  *  +  u^  (  *  +/s  o<^  ^  hrf  uCf  ^  +  ?  ij )    (1-22) 

In  Ref.  1,  E.  G.  Brown-Edwards  calculated  the  aerodynamic 
coefficients  for  a  typical  missile  configuration.   This  missile 
is  sketched  in  Fig.  2.   In  Appendix  I  Brown  Edwards*  results  are 
given  and  are  used  as  the  basis  of  the  calculation  of  the  constants 
in  Eq.  1-22.   The  body  diameter  was  taken  as  one  foot  and  the  flight 
Mach  number  as  three.   The  calculations  for  a  flight  altitude  of 


as  follows: 

•*• 

- 

- 

I 

r, 

= 

- 

z« 

*  \ 

2mV  "J^v( 

Ki 

= 

M 
1 

i 

•*, 

z 

- 

Zn 

mV^  mj£( 

& 

= 

- 

M^ 
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50,000'  are  shown  in  the  appendix.   The  results  are  tabulated  in 
Table  I.   This  table  also  contains  the  results  for  sea  level 
flight,  the  calculations  for  which  are  noir  shown. 

SECTION  2 :  THE  TRANSIENT  RESPONSE  OF  THE  MISSILE  TO  A 
STEP  DISPLACEMENT  OF  ELEVATOR  ANGLE 

Equation  1-22  relates  the  angle-of-attack  °<  to  the  elevator 
angle   ^  .   We  wish  to  find  the  variation  of   <=<   with  time  upon 
the  application  of  a  step  displacement  in  >\    ,  initial  values  all 
being  zero.   No  analytical  solution  for  this  non-linear  problem 
exists;  an  iterative  method  will  be  sought. 

Equation  1-22  is,  in  operational  form, 

(D2  +  2  5^  D  +  uC,  )  ^  +/3«<5  =  Kx  ul^d  +  z  D)  ^       t2"1) 
This  may  be  written: 

T(D)  K  +  p>   ^   =  U(D)^  (2-2) 

where 

T(D)  =   D2  +  2  T,  ">~t  D+  u)^(  (2-3) 

U(D)  =   Kx  uJ^  (l  +  Z   D)  (2-4) 

Now  an  approximation  to  the  solution  of  Eq.  2-2  is  the 
solution  o(x  (t)  of 

T(D)  *,   =  U(D)  ^  (2-5) 

This  equation  is  linear,  and  could  be  solved  analytically,  however 
we  will  use  an  approximate  method.   In  Ref .  1,  where  an  attempt  is 
made  to  solve  the  same  non-linear  problem,  Eq.  2-5  .is  solved 
analytically.   When  this  is  done  the  ensuing  steps  become  quite 
tedious.   We  hope  to  find  an  advantage  in  taking  another  course. 
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We  find   <^,  (t)  in  the  form  of  a  Fourier  series  using  the 
approximate  method  given  by  Wass  and  Haymen  in  Ref .  2  and  described 
in  Appendix  II  of  this  thesis.   The  method  requires  that  the 
frequency  response  of  Eq.  2-5  be  determined.   The  reader  who  is 
unfamiliar  with  frequency  response  technique  will  find  it  explained 
in  Section  5»  in  connection  with  the  design  of  the  control  system. 
The  Wass,  Hayman  method  gives   c<  (t), 
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approximately,   in  the  form, 

*,<*>  ■  hi  [*%  +  4  Z«  sin  (<n^  *  +  B"}"1        (2_6) 

>>v  =  1,  3,  •  •  •  »11. 

where  u)  is  a  known  constant,  |»^|  is  the  magnitude  of  the  step 
input,  and  the  A's  and  B's  are  known  constants  -  gains  and  phase 
shifts  of  the  frequency  response.   If  positive  damping  exists, 
this  response  is  a  damped  oscillation,  with  a  steady-state  value  of 


(*,)B  =  WK  (2-7) 

The  true  steady-state  value  (  <K  )  may  be  calculated  from  Eq.  1-22 
by  setting  all  the  terms  involving  derivatives  equal  to  zero: 

The  solution  for  (  <K  )s  is  easy,  though  not  analytic. 
Now  if  we  take,  as  a  new  approximation  to  o< (t), 


*2(t) 

= 

(c<)s 

o<±(t) 

= 

A 

0 

[*• 

2-I- 

then 

its 

steady- 

•state  value 

is 

(   °< 

K        2 

:>„ 

A 

0 

■      Ao    ■ 

■    C0B 

(2-9) 


-Since  °v,(t)  has  the  true  steady-state  value,  there  is  some 
reason  for  thinking  it  is  an  improved  approximation  for  o<(t). 
Some  further  justification  is  presented  in  Fig.  3.   The  curves  in 
this  figure  represent  the  solutions,  by  digital  computer,  of  an 
equation  essentially  the  same  as  (l-22).   The  equation  is  discussed 
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in  detail  in  Section  3.   The  solutions  for  the  linear  case  and  for 
several  non-zero  values  of  the  coefficient  in  the  cubic  term  are 
plotted.   In  relation  to  the  present  example,  the  linear  solution 
corresponds  to  °^-,(t),  and  the  non-linear  solutions  to  o<   (t). 

In  Fig.  4,  each  of  the  non-linear  solutions  has  been  scaled  up  to 
give  steady-state  values  equal  to  ( °^,  )  -  an  inversion  of  the 
process  represented  by  Eq.  2-9.   It  may  be  seen  that  this  scaling 
process  apparently  brings  the  non-linear  solutions  into  closer 
agreement  with  the  linear  one  in  a  general  sense,  not  just  the  steady- 
state  values.   Thus  °^At),   a  scaled-down  version  of  °*-,\.t)t   will 
be  taken  as  an  improved  approximation  to  ©<(t).   For  brevity  let 

R  =  v  ;s  ^2-10) 

h«Ao 
so  that 

«*2(t)-     R^W  (2-11) 

For  later  convenience  let 

**!    (t)  =    ^o  +  £    V^nsin  ("lULt  +  B J  (2-12) 

i.e.  the  response  of  the  linear  system  to  a  unit  step  input. 
Thus, 

°Vt;  =    hi  ^  (t)  (2-13) 


and 


<*2(t)  =    Rh|  oi±   (t)  (2-14) 


Now  if  /3  is  sufficiently  small,  a  closer  approximation 
to  c*(t)  is  the  solution  oi(\)   of 

T(ft)o<3  =  u(D)rj   -  /S^ 

=  U(D)7  -  ft    r5  |^|  3  o<  3  (2_15 


Or,   if      °^^^(t)  is  the  solution  of 
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T(D)*V=     -£<*2=-/Sr3    hi3     *1 


then  by  the  principle  of  linear  superposition, 


(2-16) 


*5(t)  =   ^x(t)  +  o<5#(t)  (2-17) 

To  solve  (2-16),  a  method  will  be  developed  for  expressing 


oi  ,    as  a  Fourier  series. 


A~  sin  (*vlO  t       )" 


1 


•n  =  l,  3,  . . . ,  ii. 

The  expanded  cubic  consists  of  the  cube  of  each  of  the  original 
terms,  plus  three  times  the  product  of  the  square  of  each  term  and 
each  of  the  other  terms,  plus  six  times  the  product  of  all 
combinations  of  three  unlike  terms.   Showing  a  few  terms  of  the 
expansion  and  introducing  a  notation  for  their  coefficients  we 
have: 

1»      "  a000  ko     +  a001  Ao2  Al  Sin  (t°Ft  +  V  +  ••' 

+  am  ki  sin3  (%* +  v 

+    aU3  Ax2  A3  sin2   (u>pt  +  B±)  sin  (3**>pt  +  B^)  +  ... 


+     a.  ..    A. A  .A,    sin  (iu)_,t  +  B.)  sin  (ju),,t  +  B.)  sin  (kui_t  +  B,  ) 
xjk     l  j  k  v       F  i'  Kt}     F  j'  x       F  k' 


(2-18) 


Each  term  is  identified  by  a  subscript  i,  j,  k.  The  coefficients 
a.  .,  take  account  of  the  kind  of  term  (e.g.  the  multiplier  6  for  a 
term  which  is  the  product  of  three  unlike  terms)  and  the  factors 
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l/2,  2/tt  »  and  l/n.   The  coefficients  a.  ..  are  tabulated  in  Table  II. 
Any  product  of  sine  terms  can  be  expressed  in  the  form, 


E 


o 


+  E,  sin  (iO„t  +  P. )  +  E0  sin  (2  iowt  +  F_) 
1       F     12        F     2 


+   ...   +  E^  sin  ( a  «Opt  +  F^  )  +  . . . 
For  example,  in  the  1,3,5  term  we  have  the  identity: 
sin  (u)pt  +  B^  sin  (3«Opt  +  B^)  sin  (5«Opt  +  B^) 

a  -  i  sin  (u)pt  -B1-B5+B5)  +  i  sin  (3^pt  +  Bj-B  +B5) 

+  i  sin  (7u>pt  -B1+B3+B5)  -  i  sin  (9u>pt  +B]+B3+B5) 

In  this  case 

IP  _         _   J.  TT—        J-  T?  _        J-  TT_        L 

JJj^        -        -     4  ,  JJi,        -       "4"  ,  £i^       -        4      ,  Jig        -        -4- 

E0  =  E2  =  E4=E5  =  B6  =  E8  =  E10=     ...     =     0 

r1-   -Bi-b3  +  b5 

F~  =     B,    —  B^  +  Br- 
3  13        5 

P„  =     -  B,   +  B^  +  Bc 
7  13         5 

F9  =    B!  +  B3  +  B5 

Using  this  notation  a  large  number  of  these  identities  may  be 
tabulated  concisely,   and  this  has  been  done,   for  the  combinations 
necessary  here,   in  Table  III. 

We  adopt  the  notation*: 
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sin  (iu^t  +  Bj.)  sin  (j  u>?t  +  B^)  sin  (kw>pt  +  B^ 

=  Eo(ijk)  +  El(ijk)  Sin  ("V  +  Fl(ijk)}  +  •'•  <*■«> 

Also,  we  will  write,  for  later  convenience, 

sin  (iU)pt  +  Bi)  sin  (jiopt  +  B^) 

=  E  /  . . \  +  E, /  . . \  sin  ( to  t  +  F, /  .  .  x )  +  . . .  (2-20) 

o(oijj    l(oij)       F    l(oijr 

and 

sin  (iu>pt  +  B±) 

=  E  /   •  \  +  E,  /   .  \  sin  (toj;  +  F,  /   .J  +  ...  (2-2l) 

0(001;    l(ooiJ     v  F     1(001;'  x    ' 

The  latter  expression  is  a  trivial  use  of  such  a  notation  since  the 
right  hand  side  reduces  to  sin  (i*o >  t  +  B.)»  but  its  use  helps 

£  1 

provide  a  systematic  tabulation. 

Using  the  identities  represented  by  Eqs.  2-19,  2-20,  and  2-21 

in  Eq.  2-18  gives  °^ ,    as  a  series  of  constants  and  simple  sine 

1  u 

terms.   Since  sine  terms  of  like  frequency  can  be  combined  into  a 
single  sine  term  of  that  frequency,  we  have: 

*1  5  =  bo  +  bl  sin  ("V  +  Cl)  +  b2  sin  ^"f*  +  °2)  (2*22) 

+  ...  +  b   sin  (  /vl  to  t  +  C  „  )  +  . . . 
<»t         F 

where 
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3        2 
o    000  o     001  o   1 


+  a. 


i.  ..  E  /.  ..  \  A.  A.  A.  +  ... 

ljk  o(ijk;  l  j  Tc 

bl  sin  H*  +  Cl}  "  a001  E1(001)  Ao2  Al  Sin  ^P*  +  F1(001)} 
+  a003  El(003)  Ao2  A3  Sin  (uV  +  Fl(003)} 

+  aijk  El(ijk)  Ai  Aj  \  Sin  (l°F  t+  Fl(ijk)>  +  '•• 

b   sin  (n  vA>T,t  +  C  )  =  ... 
n         F     n 

+  aljk  En(i*)  Ai  Aj  ^  sin  <n  -V*  +  *n(uk))  +--- 

It  is  possible  to  shorten  the  notation  conveniently  by 
letting 


nUjkj     ljk   n(ijk)  v   v/ 

The  values  of  c  / .  .,  \  may  be  determined  from  the  tabulations  of 
aijk  and  En(ijk)*   Accordingly,  values  of  cn(-Hk)  (or  simply  c) 
are  tabulated  in  Table  IV.   The  coefficients  are  all  listed  as 
positive.   The  coefficients  for  which  the  values  of  E  were  negative 
were  changed  in  sign  by  adding  180  to  the  corresponding  angles  F, 
which  are  also  tabulated.   It  will  be  seen  that  like  harmonics 
have  been  grouped  together. 

Table  IV  does  not  include  coefficients  for  all  possible  terms. 
The  elimination  of  less  important  terms  was  based  on  the  effect  of 
each  term  on  °^(t)  (the  approximation  towards  which  we  are  working) 
in  the  numerical  example  which  follows.   It  had  been  hoped,  at  the 
outset  of  this  investigation,  that  the  number  of  terms  needed  would 
be  small.   But  it  was  found  necessary  to  retain  a  large  number  of 
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terms,  up  to  and  including  9th  harmonic  terms.   The  relative 
importance  of  the  smallest  terms  retained  can  best  be  seen  in  the 
example  of  this  section. 

The  fact  that  the  requirements  of  a  particular  numerical 
example  must  be  injected  at  this  point  is  a  deficiency  in  the  method, 
since  the  choice  of  retained  terms  is  not  necessarily  the  best  choice 
for  other  examples  to  which  this  method  will  be  applied.   However, 
to  follow  a  more  general  course  would  be  impractical. 

To  obtain  values  for  b   and  C   for  a  particular  problem, 
the  products  c  / .  .,  \  A.  A.  A,   and  the  angles  F  / .  .,  \  are 
calculated  using  the  values  of  gain  A  and  phase  shift  B  from  the 
frequency  response  curves,  then  like  harmonics  are  combined.   The 
combination  of  several  sine  terms  in  the  same  frequency  may  be  done 
most  easily  by  graphical  vector  addition.   This  is  described  in 
Appendix  III. 

The  Fourier  series  for  °^,    may  be  computed  with  minimum 

1  v 
effort  by  using  the  form  of  Table  IV.   All  the  necessary  results 

of  the  preceding  development  are  summarised  there,  and  the  required 

computations  are  indicated  by  its  layout.   To  use  the  table, 

proceed  as  follows: 

1.   Enter  the  values  of  gain  A  and  phase  shift  B  of  the 
appropriate  frequency  response.   Once  this  is  done,  no 
further  reference  need  be  made  to  any  information  other 
than  that  contained  in  Table  IV. 


2.  Calculate  the  products  A.  A.  A,  . 

3.  Calculate  the  products  c  .  .,  \  A.  A.  A,.   (in  the 

table   cAAA  is  written  for  c  / .  .,  \  A.  A.  A  ) 

n(ijk;   l  j  k' 
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4.  Calculate  the  angles  F  /.  .,  \. 

5.  Combine  like  harmonics  (e.g.  by  the  method  of 
Appendix  III),  and  enter  the  results,  bn  and  Cn  , 
in  the  space  provided. 

Step  5  above  gives  the  Fourier  series  for  °^ 1      .   Multiplying 
each  of  the  coefficients  of  the  series  by  the  constant  -  R  |*l| 
gives  the  Fourier  series  for  -  IT  |^|   c/   .   By  Eq.  2-16, 
o<  ^*(t)  may  be  found  by  applying  the  frequency  response  P/t(x«o). 

The  next  approximation  to  °<(t),  if  it  is  desired  to  go 
further,  is  °C,(t),  the  solution  of 

T(B)o^4  =U(D)  1   -/3*53  (2-24) 

The  determination  of  a  Fourier  series  for  °^_   may  be  obviated 

•z    <s 

as  follows.   Adding  and  subtracting  p  °*    on  the  right  hand 
side  of  (2-24): 

T(D)o<4  =  u(D)ij  -^^23  +  /3  (<*/  -  <*/)      (2-25) 

Let     °C,#(t)     be  the  solution  of 

T(D)  o<^    =    £  (<X23  _    o<35)  (2-26) 

Then  by  superposition  of  the  solutions  of  Eqs.  2-15  and  2-26, 


cxf4(t)  =   ©<5(t)  +  <*4»(t)  (2-27) 

The  Fourier  series  for  (°^?  -  °C  )  need  not  be  determined  with 
as  great  an  accuracy  as  would  be  required  for  that  of  °C,  .   To 
solve  Eq.  2-26,  °<  and  °C,  are  each  evaluated  and  cubed,  point 
by  point.   The  difference  is  then  plotted,  and  a  numerical  harmonic 
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analysis  is  performed.   The  frequency  response  ^/T(ioo)  is  applied 

to  the  resulting  courier  series,  giving  ^.^       Instead  of 

4* 

evaluating  ^C,(t)  and  cubing  each  point,  the  series  which  has  been    -j-  j 

Shu  Act****.  <~utk*l  <~  p^^oU/^t-j  '*/  °<^  Acvt,  t*+vr  jiA/ixl^a^2^-, 

obtained  for  °v>  could  oe  evaluated, A be cause  cubing  a  few  numbers 

is  much  easier  than  evaluating  a  Fourier  series.   But  the  latter 
method  has  a  certain  superiority  in  that  if  the  series  obtained  for 
°(.  „  is  not  a  good  representation  of  the  true  cube  of  °<  ,  due 
either  to  a  mistake  or  to  a  deficiency  in  the  cubing  method,  then 
the  error  will  be  compensated  only  if  the  latter  method  is  used. 
A  deficiency  in  the  cubing  method  would  probably  exist  if  the  linear 
system  frequency  response  was  very  different  from  that  used  in 
formulating  the  cubing  method. 

Improved  approximations,  °C,  <?C,  etc.,  would  be  found  in  the 
same  way  as  is  °(.. 

A  method  is  shown  in  Section  6,  in  connection  with  the  control 
system  response,  for  simplifying  the  calculations  when  the  responses 
to  more  than  one  input  magnitude  are  required.   An  analogous 
method  can  and  should  be  applied  here  if  more  than  one  input 
magnitude  is  used. 

We  now  apply  this  method  for  the  case  of  missile  flight  at  sea 
level  and  an  input  magnitude  of  .25  radian.   This  case  was  chosen 
because  its  analog  computer  solution  was  available.   The  constants 
m  Eqs.^2-3,  and  2-4  are  those  tabulated  for  sea  level  in  Table  I, 
with  one  exception:  *Z     is  taken  as  zero.   The  effect  of  that 
term  on  the  transient  response  was  found  to  be  negligible,  having 
an  influence  on  the  frequency  response  only  at  very  high  frequencies, 
and  it  was  omitted  from  the  analog  computer  setup.   It  was 
retained,  in  the  development  just  shown,  for  generality. 
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The  frequency  response  of  Eq.  2-5, 

U(iu>)     =  X  1 


(iu>f  +  2|5X  u;^  (ito)    +  uuA 


w£\*  +     ^ 


Kl (2-28) 


(*) 


is  plotted  in  Fig.  5.   Gains  are  plotted  in  decibels,  and  phase 
shifts  in  degrees.   Following  the  Wass,  Hayman  method  for  determining 
linear  transient  response,  we  take  UX,  to  be  one  fifth  of  the 
frequency  at  peak  gain.   The  peak  occurs  at 


CO  =   LO^ 


So  that 


u)   =      2  tO 


The  gains  A   and  phase  shifts  B   of  Eq.  2-28,  at  frequencies 
nu)p  where  n  =  0,  1,  2,  ...,  13,  are  tabulated  in  Columns  2  and 
3  of  Table  V,  with  values  of  n  in  Column  1.   The  values  of  gain 
are  absolute,  not  decibels.   This  procedure  is  necessary  to  make 
the  computations  systematic  (this  differs  slightly  from  the  Wass, 
Hayman  procedure). 
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The  coefficients  of  the  sine  terms  in  the  Fourier  series'  for 
a  square  wave,  magnitude  .25,  are  tabulated  in  Column  4  of  Table  V. 
Each  coefficient  is  identified  with  the  harmonic  to  which  it  applies 
by  the  value  of  n  in  the  same  row,  e.g.  the  coefficient  of  the 
sin  7  ^t  term,  or  seventh  harmonic,  is  .0228.   This  procedure  will 
be  followed  throughout  the  thesis,  including  the  use  of  a  second 
column  for  phase  angle  when  required,  and  will  be  referred  to  as 
"tabulating  the  Fourier  series." 

Having  tabulated  the  Fourier  series  for  the  input  square  wave, 
and  the  frequency  response  of  the  linear  system,  we  find  the 
coefficients  of  the  Fourier  series  for  the  output,  °^, (t),  by 
multiplying  each  input  coefficient  by  the  corresponding'  gain. 
The  pnase  angles  of  the  output  series  are,  in  this  case,  the  phase 
shifts  of  the  frequency  response.   The  Fourier  series  for  °*, (t) 
is  tabulated  in  Columns  5  and  6.   From  this  tabulation  we  see,  for 
instance,  that  the  fifth  harmonic  term  of  the  Fourier  series  for 
o<1(t)  is  .0775  sin  (5**>Ft  -  88.1°). 

The  evaluation  of  <K,(t),  for  39  values  of  time  up  to  95 
percent  of  the  half-period  of  the  square  wave,  is  tabulated  in 
Table  VI.   The  evaluation  was  made,  as  were  all  of  those  in  this 
report,  using  the  Wass,  Hayman  method.  °<  ,(t)  is  plotted  in 
Fig.  6  together  with  the  corresponding  solution  obtained  by 
Brown-Edwards  on  a  Short  Analog  Computer  (see  Ref.  l).   The 
oscillation  is  seen  to  be  of  high  amplitude,  due  to  the  low  damping. 
It  is  perhaps  surprising  that  the  approximation  is  as  good  as  it  is, 
because  the  Wass,  Hayman  method  assumes  that  the  oscillation  will  be 
nearly  completely  damped  out  by  the  end  of  the  half  cycle  of  the 
square  wave,  and  it  clearly  is  not. 

Solving  Eq.  2-8  for   |w  |  =  .25  we  find 
(«<)s   =  .195 
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By  Eq.  2-10, 

R  =     t.25J(%S3>    "  -917 

By  Eq.  2-11,  the  evaluation  of  °^2^)  is  found-  from  the  evaluation 
of  °^-,(t).   This  is  tabulated  in  Table  VI,  and  plotted,  with  the 
analog  computer  solution  of  the  non-linear  equation,  in  Fig.  7. 

©< iy  is  calculated  by  the  method  of  Table  IV.   This  calculation 
is  shown  in  Table  VII,  and  the  resulting  Fourier  series  has  been 
transferred  to  Columns  7  and  8  of  Table  V.   Multiplying  the 
magnitudes  by 

R3(»j|5  =  (.917)5  (.25)3  =  .01206 

and  adding  180  to  the  phase  angles  gives  the  series  for 
-  R  |»]  |   <*\   ,  or  -  °<   ,  and  this  is  tabulated  in  Columns  9 
and  10.   The  frequency  response  /^/T(iu>),  which  in  this  case  only 
differs  from  U(i**))/T(i  u>)  by  a  constant  multiplier,  is  tabulated 
in  Columns  11  and  12.   By  Eq.  2-16,  °<  _  is  found  by  applying  the 
frequency  response  r/^{x^)     to  the  Fourier  series  for  -  <*t~   , 
and  the  resulting  Fourier  series  is  tabulated  in  Columns  13  and  14. 
This  is  evaluated  for  27  values  of  time  in  Table  VI,  and,  by 
Eq.  2-17,  combined  with  ©^(t)  to  give  °C(t).    °<5(t)  is 
plotted  in  Fig.  7. 

In  Fig.  7  we  see  the  relationship  between  the  "reduced"  linear 
solution  °^p(t),  the  result  of  the  first  iteration  ^^(t),  and  the 
computer  solution  ©^(t).   Superficially  ©^„  doesn't  look  much 
better  than  °^p,  particularly  with  respect  to  the  amplitudes  of 
the  peaks  and  troughs.   The  times  at  which  the  overshoots  occur 
have  been  brought  into  line  however,  and  it  will  be  shown  in  Section  3 
that  the  next  iteration  gives  rapid  convergence.   The  solution  to 
the  present  example  was  not  taken  further  because  the  digital  computer 


25 


solution  of  a  similar  equation  became  available  to  the  author,  and 
the  analog  computer  solution  was  thought  not  to  have  sufficient 
accuracy  to  justify  its  use  as  a  standard  with  which  to  compare 
more  refined  approximations. 

SECTION  3:  COMPARISON  01?'  RESULTS  WITH  SOLUTION  OBTAINED  WITH 
A  DIGITAL  COMPUTER 

Mr.  P.  A.  T.  Christopher  obtained  solutions,  from  a  Perranti 
Mercury  digital  computer,  of  the  equation: 

2  3 

D,  x  +  a,Dx  +  a  x  +  bx   =  Q  (3-1 ) 

where 


al 

= 

•  6, 

a 

0 

= 

1, 

Q  is 

a 

step  function, 

Initial 

values 

are  zero. 

These  solutions  were  obtained  for  various  values  of  b  and  magnitudes 
of  Q.   The  accuracy  of  the  solutions  was  such  that,  for  present 
purposes,  they  may  be  considered  exact. 

The  output  x(t),  of  Eq.  3-1,  is  dependent  on  four  parameters: 
a,,  a  ,  b,  and  1  Q J  ,  the  magnitude  of  the  step  input.   It  is  shown 
in  Appendix  IV  that  by  linear  transformations  of  the  variables  x 
and  t  ,  the  output  can  be  made  a  function  of  only  two  parameters: 
the  damping  ratio  £*  »  an^-   the  non-dimensional  input  magnitude  |Q  J  . 
Thus  the  solution  of  Eq.  1-22,  which  is  of  the  same  form  when 
"Z   =  0,  could  be  found  from  the  solution  to  Eq.  3-1,  if  V  and 
|  Q  |   for  each  of  the  equations  were  identical.   Now  the  solutions 
to  Eq.  3-1  include  a  range  of  values  of  |Q0|  ,  but  only  one  value 
of  $   ,  i.e.   .3,  while  the  value  of  ^  for  the  missile  at  sea  level 
is  .175.   Accepting  this  deficiency,  we  will  apply  the  approximation 
method  to  Eq.  3-1  for  the  computer  value  of  |Q  1   nearest  that  of 
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the  case  considered  in  Section  2.   The  results  will  not  be 
transformed  to  the  amplitude  and  time  scales  of  the  missile  problem, 
because  having  already  seen  the  approximate  missile  solution  in 
Section  2,  we  are  now  only  interested  in  the  convergence  of  the 
approximations  -  which  does  not  depend  upon  scale. 

The  equivalent  symbols  in  Eqs.  3-1  and  1-22  are  tabulated 
below  (7=0): 

Eq.  3-1  Eq.  1-22 

2 


o 
b 


~1 


Thus  the  non-dimensional  input  in  terms  relating  to  Eq.  1-22  is: 


^5    =   ^1^-,  1)J£    =    Kl^  1 


^o        V2      ■L 


5/2 

(U)2  )  *2 


2~W; 


u) 


'W 


1 


The  value  of  |Q  |   corresponding  to  the  case  considered  in  Section  2, 
where  the  constants  are  those  given  in  Table  I  for  sea  level  flight 
and    \v\\  =   .25,  is: 


\Q  \    =   (-855)r(l546)(.25)    =  r^ 

°  (24.4) 

The  computer  solutions  nearest  this  are: 

(l)   b    =  .5 

Ul   =  .4   __ 
\%\     =      *±£Z         =  .283 
(1)3/2 
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(2)   b    =1 

tal  =  .4 

|Q  1  =    .4  fl      =     .400 

We  choose  the  higher  value,  case  (2),  and  proceed  with  the  application 
of  the  approximation  method  to  this  case.   In  order  to  use  the 
equations  of  Section  2  unaltered,  we  put  Eq.  3-1  in  the  form  of 
Eq.  2-2  by  letting 

ot      =  x 

/3   =  b  =  1 

T(D)  =  D2  +  a,D  +  a   =  D2  +  .6D  +  1 
1     o 

U(D)  =  a   =  1 
x      o 

o 
llj  =  kl  =  .4 

The  gain  of  the  frequency  response 


U(icU)       

T(itL>)       1  -  uj^  +  .6  uJi 


2 


(3-2) 


has  a  peak  at  tO  =  1,  so  we  take  oJu  =  .2.   Values  of  gain  and 
phase  shift  of  (3-2)  are  tabulated  in  Columns  2  and  3  of  Table  VIII, 
for  the  required  values  of  a *-<"}.!•   The  Fourier  series  for  a  square 
wave,  magnitude  .4,  is  tabulated  in  Column  4.   By  Eq.  2-5  we 
find  the  Fourier  series  for  °^(t),  and  this  is  tabulated  in  Columns 
5  and  6.   The  evaluation  of  °^, it)  for  22  values  of  time  is 
tabulated  in  Table  IX,  and  is  plotted  in  Fig.  8  with  the  corresponding 
computer  solution.   The  calculated  solution  <X,  appears  to  be  a 
good  approximation  to  the  true  output  of  the  linear  system. 
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The  steady-state  output  is  found  from  Eq.  2-2  by  setting  the 
derivative  terms  equal  to  zero: 


For  \\r\  |  =  .4  this  equation  yields 


(o<)s  =  .355 


By  Eq.  2-10, 


R  =  /  ffi?.,  \  =  .888 


The  calculation  of  °<         is  performed  in  Table  X.   The 


The  evaluation  of  °C(t)  is  found  from  the  evaluation  of  °*-,(t) 
according  to  Eq.  2-11.   This  is  tabulated  in  Table  IX,  and  is 
plotted  in  Fig.  9. 

lv 

resulting  Fourier  series  is  entered  in  Columns  7  and  8  of  Table  VIII. 

Multiplying  the  magnitudes  by 

R3  \y\\3     =     (.888)3  (.4)3  =  .0447 

and  adding  180  to  the  phase  angles  gives  the  series  for 

■z     -z       "Z  -x 

-S  1*1  I   °^n    or  -  o<^  ,  and  this  is  tabulated  in  Columns  9 
and  10.   The  frequency  response  fi/l(±uj)     is  the  same  as 
U(iu))/T(iu))  in  this  case  because  (i  =  u(i**>J  =  1.   This  has 
already  been  tabulated  in  Columns  2  and  3.   Thus  °t-z*   is  found 
according  to  Eq.  2-16,  tabulated  in  Columns  11  and  12,  and 
evaluated  in  Table  IX.  °^~   is  obtained  using  Eq.  2-17,  and  plotted 
in  Fig.  9. 

The  next  approximation  requires  a  numerical  harmonic  analysis. 

■2  'Z 

In  Table  IX,  o(      and  °<.,  are  cubed,  and  the  difference  (  <*  i     -  c<;) 
2       3  3     3  5 

is  found  for  each  of  the  22  values  of  time.   (o<   -  c<  )  is 

plotted  in  Fig.  10. 
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Wow  harmonic  analysis  requires  that  the  function  being  analysed 
be  known  throughout  one  period  of  the  fundamental  frequency  tO  . 

r 

So  far  only  the  first  half-cycles  have  been  evaluated.   The 

i^i  Xht*  <&c+*<l  ktvU,-  &$&£<-. 
evaluations  of  °^0  and  ©^  are  given,  for  22  values  of  time,  in 

Table  XI.   In  this  tabulation,  the  time  t-,  is  the  time  measured 

3      3 
from  the  start  of  the  second  half -cycle.   ©<   ,  o(      f  and 

3        3  \ 
(  °t         -  °^v  )  are  also  tabulated.   We  see  that  the  values  of 

2       3  )  in  the  second  half-cycle  are  negligible  compared 

with  those  in  the  first  half-cycle,  and  set  them  all  equal  to  zero 

in  the  numerical  harmonic  analysis.   This  result  is  not  surprising, 

because  the  oscillation  in  the  second  half-cycle  takes  place  about 

the  steady-state  value  of  zero.   Therefore  the  effect  of  the 

non-linearity  is  small,  and  the  oscillation  of  the  non-linear  system 

is  nearly  the  same  as  that  of  the  linear  one.   This  can  probably  be 

,3     .,3 

c 

evaluated  for  the  second  half -cycle. 


taken  as  a  general  result,  so  that  ( d        -      o<      )  need  not  be 


3       3  \ 
A  24  ordinate  harmonic  analysis  of  (  o^p   -  o<   ),  performed 

using  the  layout  recommended  by  Wylie  in  Ref.  3,  is  shown  in 

Table  XII.   The  positions  at  which  the  ordinates  were  taken  is 

shown  by  the  vertical  lines  in  Pig.  10.   They  were  found  by 

dividing  the  first  half-cycle  into  12  equal  time  intervals. 

The  ordinates  are  labeled  y,,  y„,  etc.,  in  accordance  with  Wylie* s 

notation.   The  ordinates  were  all  multiplied  by  1000  to  make  the 

harmonic  analysis  more  manageable.   The  result  of  the  analysis  is 

a  series  in  both  sine  and  cosine  terms  which  includes  the  12th 

harmonic.   Retaining  only  up  to  the  9th  harmonic,  the  series  was 

put  into  the  form  containing  only  sine  terms,  and  this  is  tabulated 

in  Columns  13  and  14  of  Table  VIII.   °< .#  is  obtained  by  Eq.  2-26, 

and  this  is  tabulated  in  Columns  15  and  16.   It  is  evaluated  in 

Table  IX  and,  in  accordance  with  Eq.  2-27,  added  to  <*( „     to  give 

o< ..   This  is  plotted  in  Fig.  9  with  the  other  solutions.   We 
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see  that  the  convergence  is  virtually  complete  in  the  region  of  the 
first  peak,  and  probably  good  enough  for  most  practical  purposes  in 
the  entire  range  of  the  computer  solution. 

This  completes  the  demonstration  of  the  method  developed  in 
Section  2.   It  would  be  desirable  to  find  out  how  many  iterations 
are  necessary  for  larger  values  of  input  magnitude,  and  for  what 
values  of  input  magnitude  the  first  iteration  is  sufficient,  but 
this  could  not  be  done  by  the  present  author  in  the  time  available. 

SECTION  4:  EQUATIONS  OF  COMPONENTS  Off  THE  MISSILE  CONTROL  SYSTEM 

A.   Incidence  Vane 

One  of  the  components  of  the  control  system,  which  will  be 
designed  in  Section  5,  is  an  incidence  measuring  vane  at  the  nose 
of  the  missile.   The  design  of  a  conical  type  vane  suitable  for 
this  use  is  shown  in  Appendix  V.   The  vane  will  be  pivoted  at  its 
e.g.  position. 

The  relationship  between  the  angle  of  attack  of  the  missile 
and  the  output  of  the  vane  must  be  determined.   Fig.  11  shows  the 
angles  and  sign  conventions  used  in  the  analysis.   The  angle  )f  is 
measured  by  a  transducer  in  the  vane  and  is  the  output  of  the  vane. 

It  will  be  seen  from  the  figure  that  no  distinction  is  made 
between  the  velocity  of  the  e.g.  of  the  missile  and  the  velocity  of 
the  vane.   This  is  not  strictly  correct,  of  course,  because  the 
pitching  velocity  of  the  missile  q  creates  a  component  of  vane 
velocity  perpendicular  to  the  missile  axis  with  magnitude  ly.  q. 
It  is  easily  shown,  however,  that  for  the  configuration  and  flight 
conditions  under  consideration,  this  component  of  vane  velocity  has 
negligible  effect  on  the  vane  velocity.   For  the  example  being 
considered,  under  steady-state  conditions  at  an  angle  of  attack  of 
40  ,  the  component  of  vane  velocity  due  to  pitching  of  the  missile 
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is  5.3  ft/ sec,  and  the  difference  between  vane  incidence  and  missile 


incidence  is  .12  . 


The  equation  of  angular  motion  of  the  vane  about  its  e.g.  is: 

\  =  h    ^V  (4-1) 

The  aerodynamic  moment  on  the  vane  is  taken  as 

"V  =  ci  *v  +  c2  *  (4"2) 


where 


ci  - 


a  x 


v 


9* 


and  C?  K    is  the  internal  damping  moment  arising  from  the 
relative  angular  motion  between  the  vane  and  the  missile.   It  is 
assumed  that  a  device  within  the  vane  exists  to  provide  this  damping. 

As  in  the  case  of  the  missile  (see  remarks  following  Eq.  1-3) 
the  aerodynamic  damping  moment,  that  due  to  °^v»  is  omitted  on 
the  assumption  that  it  is  small,  but  for  the  real  reason  that  a 
value  for  it  is  not  available. 

Combining  Eqs.  4-1  and  4-2  gives: 

C1^v  +  c2  V   =  Iy  ©^  (4-3) 

rte  wish  to  reduce  Eq.  4-3  to  a  linear  equation  in  the  two 
variables  o<  and  %   .   From  the  figure, 

^v  =  c*-  K  (4-4) 


and 


ev  =    eM    -  y  (4-5) 
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Substituting  (4-4)  and  (4-5)  in  (4-3)  gives 

cx  (•<-  %)  +  c2  }   =  iy  (eM-  tf  )  (4_6) 

Now   ©.,  or  q  represents  the  pitching  acceleration  of  the  missile, 
and   ^   that  of  the  vane  relative  to  the  missile.   It  is  reasonable 
to  assume,  and  it  will  be  evident  later,  that  the  frequency  of 
oscillation  of  the  vane  is  much  higher  than  that  of  the  missile. 
This  being  the  case,  and  providing  that  the  magnitude  of  the  missile 
oscillation  is  not  excessive,  then  0     may  be  assumed  negligible 
in  comparison  with   V  .   Accordingly,  Eq.  4-6  becomes 

c1  (-<-  y)  +  c2  y   =  - 1  % 


or 


v  +  -I  i  -  ^  i  =  -  -±  *  (4-7) 

v        v  v 

The  constants  u^      and   5"0  are  now  defined  as  follows: 

2        d 


2 

u)S     =  -  c 

2 


"Vif 


^2  =   C2/2  ^ 


uJ 


^v 


Eq.  4-7  then  becomes 

V  +  2  f  U3    Y  +  u)^.   V  =  tU*f  o<       (4-8) 

^2  2  2 

The  calculation  of  the  constants  will  be  found  in  Appendix  VI. 
They  are  tabulated  in  Table  I. 
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B.   Actuator 

The  equation  for  a  typical  elevator  actuator  and  the  constants 
of  that  equation  were  obtained  from  D.  Morton,  a  student  at  the 
College  of  Aeronautics  who  had  investigated  the  subject.   The 
equation  is: 

\      +     2  ^3  "V  h    +  ^     \  "  ^     1  D   (4"9) 


j 


where   *i  ^  is  the  demanded  elevator  angle, 

V)   is  the  elevator  angle. 

The  constants  uJ^,   and   $ „  are  given  in  Table  I. 
3        ^ 


SECTION  5:  DESIGN  OF  THE  CLOSED-LOOP  CONTROL  SYSTEM 

A.   The  Linear  System  in  Transfer  Function  Form 

The  control  system  arrangement  to  be  used  is  shown  in  Fig.  12 
in  block  diagram  form.   In  this  diagram: 

o<    is  the  demanded  missile  angle  of  attack,  the  control 

system  input, 
©<R   is  the  reference  input,  equal  to  K,  times  ^  where 

the  constant  IC,  is  the  gain  of  the  linear  amplifier, 
o<   is  the  missile  incidence, 
•jj  is  the  demanded  elevator  angle, 
H   is  the  elevator  angle, 
)(   is  the  output  of  the  incidence  measuring  vane, 

and 

E=   *R-tf   =K3cKD-tf  (5-1) 

Equation  4-8  relates  °<  and   V  ,  Eq.  4-9  relates  v)     and  *\„   , 
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and  Eq.  1-22  relates  o<  and  w  .   Dropping  the  non-linearity  for 
the  present  so  that  the  design  problem  is  that  of  a  linear  system, 
Eq.  1-22  becomes 


1     1         1        x    i 


The  relation  between  E  and  \~     is,  as  yet,  unspecified. 
The  Laplace  transformation,  defined  by: 


f(s)  =    /f(t)  e_st  dt 


where  f(s)  is  the  Laplace  transform  of  f(t),  will  now  be  applied 
to  Eqs.  4-8,  4-9,  5-1,  and  5-2,  assuming  that  initial  values  of 
all  the  variables  and  their  derivatives  are  zero.   Taking  the 


Laplace  transforms  of  these  expressions  -  when  initial  values  are 

sero  - 
,2/^2 


zero  -  simply  amounts  to  substituting     s     for     d/dt     and     s       for 


d  /dt   .       These  equations  then  become,   after  rearranging, 

2 


^2  2 


•?M 


%M  *,      +     2    *,.    uJ^ 


1  3 


(5-2) 


=     V(^)     =      —       — 2- (5_3) 


2 
s  Wl^)    =      -3 z ~ o—  (5-4) 


E(*)  =        K^      <*r^)     -      tf   (*0  (5-5) 


gM    ^     y(^)  a     ^u^a  +  r^)  (5_6) 
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where   V(s),   v1j(s),  and  P(s)  are  defined  as  indicated.   The 
transform  of  the  equation  relating  E  and   l^   is  given  the 


form: 


>U)   =  ^     C(s)  (5-7) 

E(s) 


where 


C(0)   =  1 

In  servomechanism  terminology,  Y(s)  is  the  transfer  function  of 
the  vane.   Similarly  w(s),  K,  p(  S  ),  and  iv?  C(s)  are  the 
transfer  functions  of  the  actuator,  the  missile  aerodynamics,  and 
the  compensating  network-amplifier. 

Following  standard  servomechanism  technique,  the  control  system 
may  be  represented  by  the  block  diagram  in  Fig.  13. 

From  Eqs.  5-3,  5-4,  5-6,  and  5-7  (or  directly  from  Fig.  13): 


-  (s)  =  K  K  C(s)  W(s)  P(s)  V(s)  (5-8) 

E 


This  is  known  as  the  open-loop  transfer  function,  and  the  system 
comprised  of  the  four  elements  whose  transfer  functions  appear  in  it 
and  which  has  an  input  E(t)  and  an  output  %{t)     is  the 
"open-loop  system". 

Combining  Eqs.  5-3,  5-5,  and  5-8  to  eliminate  E(s)  and   V(s) 
gives  the  closed-loop  transfer  function: 

SJ  K2  K2  K^  C(s)  W(s)  P(s) 

"3j  (S)  =    I   +  %  K2  C(s)  W(s)  P(s)  7(b)  (5'9) 

C(s)  has  not  been  specified.  It  has  been  included  in  order 
to  provide  the  system  with  satisfactory  stability  characteristics. 
The  stability  of  the  system  and  the  required  form  for  C(s)  will 
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be  determined  using  frequency  response  technique.   The  choice  of 
K?  and  K_  will  first  be  discussed. 

B.   Steady-State  Accuracy 

Suppose  that  ^^  is  set  at  some  steady  value  (°< )  .   The 
rules  of  Laplace  transforms  tell  us  that  (in  this  type  of  system) 
o<   takes  on  a  steady  value  (^^  after  the  transient  motion 
has  disappeared,  and  the  ratio  (°*)c/(c*T))q  is  given  by  setting 
s  =  0  in   ^(s)/^  (s).   Accordingly,  from  Eq.  5-9: 

£V  =   f-(o)  =     ^^-  (5-io) 


since 


C(0)  =  W(0)  =  P(0)  =  v(o)  =  1 


Ideally  the  ratio  of  steady-state  values  should  be  one. 
That  is  the  steady-state  output  should  equal  the  steady-state 
demand.   We  see  from  Eq.  5-10  that  this  can  be  accomplished  by 
a  suitable  choice  of  IL»,  regardless  of  the  value  of  K,  Kp  . 
But  there  is  another  important  consideration:  if  the  missile  were 
to  operate  entirely  at  one  flight  condition,  that  is  at  a  particular 
weight,  e.g.  position,  speed,  and  altitude,  then  K-,  would  always 
have  the  same  value.   But  since  missiles  rarely  operate  under  such 
conditions,  we  may  assume  that  K,  will  not  maintain  a  single 
value  throughout  the  flight,  and  that  {c^)^^   )       will  therefore 
vary  (unless  K-,  is  made  to  automatically  compensate  it).   It  is 
evident  in  Eq.  5-10  that  the  higher  the  value  of  K,  K?  can  be 
made,  the  less  will  be  tne  variation  in  (,©<)/ (°<)   for  a  given 
change  in  K, .   K.,   is  determined  by  the  aerodynamics  of  the 
missile,  but  K?  may  be  chosen. 

Aside  from  considerations  of  amplifier  design,  there  are  two 
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factors  which  impose  an  upper  limit  on  L  K  .   One  is  that  as 
K.  Kp  is  increased,  the  system  becomes  less  stable,  and  finally 
unstable.   So  that  the  choice  of  K~  must  be  a  compromise  between 
stability  and  steady-state  accuracy.   The  other  limit  imposed  upon 
K   is  due  to  the  presence  of  spurious  inputs,  or  noise,  which 
must  not  be  amplified  excessively. 

When  K„  has  been  chosen,  IC  may  be  chosen  to  give 


=  1 


V^D;S        1  +  K  K2 

for  a  given  flight  condition. 

If  the  steady-state  equation  for  the  aerodynamics,  including 
the  non-linear  term,  is  combined  with  the  steady-state  equations 
for  the  other  system  components  to  eliminate  all  the  variables 
except  (°Os  and-  (°^t)^S*  "^e  following  equation  is  obtained: 

™*  +  -^ <<*)!  =  A5l5l  (.<  >        (5.u) 

"\  (1  +  ^  Kg  )     b      1  f  K±  K2 

We  see  that  the  larger  K,  Kp  is  made,  the  smaller  is  the 
coefficient  of  the  cubic  term.   Thus  a  high  open  loop  gain,  seen 
to  be  desirable  in  the  linear  system,  reduces  the  effective 
non-linearity  in  the  non-linear  system. 


C.   The  Frequency  Response  of  the  Linear  System 

The  frequency  response  of  a  linear  system  at  the  frequency 
uJ  is  defined  as  the  complex  number  obtained  by  substituting  i.u>  for 
s  in  the  transfer  function  relating  the  output  to  the  input.   It 
is  usually  expressed  in  polar  form,  with  the  magnitude  in  units  of 
decibels.   A  number  is  expressed  in  decibels  by  multiplying  its 
common  logarithm  by  20. 
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Now  in  a  linear  system  a  steady  sinusoidal  input  of  constant 
frequency  to  will  give  rise  to  a  steady  sinusoidal  output  with 
the  same  frequency,  after  the  transient  motion  has  disappeared. 
The  rules  of  Laplace  transforms  tell  us  that  the  ratio  of  the 
output  sinusoidal  wave  amplitude  to  the  input  sinusoidal  wave 
amplitude,  or  gain,  is  equal  to  the  magnitude  of  the  frequency 
response  of  the  system  at  frequency  u) ;  and  that  the  angle  by 
which  the  output  sinusoidal  wave  leads  the  input  sinusoidal  wave, 
or  phase  shift,  is  equal  to  the  argument  of  the  frequency  response 
at  frequency  CO  .   Furthermore,  by  virtue  of  the  properties  of 
complex  numbers,  it  can  be  shown  that  if  the  transfer  function  of 
a  system  is  the  product  of  two  or  more  individual  transfer  functions, 
then  the  phase  shift  of  the  system  is  equal  to  the  sum  of  the  phase 
shifts  of  the  individual  transfer  functions;  and,  if  the  gains  are 
expressed  in  decibels,  then  the  gain  of  the  system  is  equal  to  the 
sum  of  the  gains  of  the  individual  transfer  functions. 

The  frequency  response  curves  of  a  system,  sometimes  called 
simply  the  frequency  response,  are  graphs  of  the  system  gain  and 
phase  shift  versus  to  (or  oO   times  some  constant). 

It  will  be  seen  that  the  frequency  response  of  a  system  may  be 
found  by  substituting  iu)     for  D  in  the  differential  equation  of 
the  system,  since  the  transfer  function  is  found  by  substituting 
s  for  D  and  the  frequency  response  by  substituting  itO  for  s. 
This  procedure  is  used  elsewhere  in  this  thesis. 

As  is  shown  in  textbooks  on  servome onanisms  (see  Bibliography), 
the  frequency  response  curves  of  the  open-loop  system  determine  the 
stability  of  the  closed-loop  system,  and  they  may  be  used  as  the 
basis  of  trial-and-error  design. 

Evaluating  the  individual  frequency  response  curves  and 
combining  them  in  accordance  with  the  rules  previously  stated, 
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p(ito)  V"(i  <o)  W(i<o)  was  found  and  is  plotted  in  Fig.  14.   The 
parameter  lO/uJ—   is  used  as  the  abscissa,  and  is  plotted  on  a 
logarithmic  scale  for  convenience.   The  constants  are  those 
given  in  Table  I  for  flight  at  50,000'.   It  can  be  seen  from 
Eq.  5-8  that  P(ito)  V(iu))  W(iu>)  need  only  be  multiplied  by 
K,  K_  C(iuj)  to  give  the  open-loop  frequency  response.   It  is 
now  necessary  to  choose  a  particular  compensating  network. 

In  ordinary  servo  design  the  choice  of  the  compensator  is 
based  on  the  two  main  requirements:  l)  System  stability,  and 
2)  Independence  of  the  closed-loop  system  from  changes  in  the 
constants  of  the  components,  i.e.  high  open-loop  gain.   But  in 
the  present  example  an  additional  requirement  must  be  met:  the 
system  must  have  a  frequency  spread  amenable  to  the  Wass,  Hayman 
method  of  transient  response  calculation  (Appendix  II).   In  fact, 
it  was  at  this  point  in  the  analysis  that  the  writer  first  became 
aware  of  the  frequency  spread  requirement,  by  having  chosen  a 
compensator  that  did  not  satisfy  it.   To  demonstrate:  This  was 
a  lead-lag  network  with  the  transfer  function: 


0(s)   = 


(i  +  £-  »)(i  ♦    ijffi  o 

"V I  "Til 

(1+^2-  s)(i  +  ..125.   s  ) 


<4~ 


UJ 


Using  the  logmodang  plot,  Nichols  chart  technique  (see  Bibliography), 
the  open-loop  gain  K-,  Kp  was  set  at  17  db,  giving  a  closed-loop 
peak  gain  in  the  upper  range  of  acceptable  values  -  i.e.  a  lightly 
damped  transient.   IU  was  set  at  1.14  db  to  make  the  closed-loop 
steady-state  gain  one.   The  closed-loop  frequency  response  curves 
are  plotted  in  Fig.  15.   The  Wass,  Hayman  method  was  applied  to 
this  linear  system  to  find  the  transient  response  for  a  unit  step 
input  of  ©(^  .   The  fundamental  frequency  of  the  input  square  wave 
was  taken  as  .3  "-U  ,  one  fifth  of  the  frequency  at  the  resonance 
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peak.   The  13th  harmonic  was  included.   The  transient  response 
thus  derived  is  plotted  in  Fig.  16.   It  appears  to  settle  down  to 
a  steady-state  value  of  approximately  .O.   But  the  actual 
steady-state  value  is  one.   The  method  has  failed,  because  the 
frequency  spread  is  too  large.   We  see  from  Fig.  15  that  the 
frequency  spread,  as  defined  in  Appendix  II,  is  approximately: 


:«    =  75 


Many  other  compensating  networks  were  tried  in  an  effort  to 
obtain  a  satisfactory  compromise  between  the  three  stated 
requirements.   In  fact,  none  was  found,  and  it  appears  to  the 
writer  that  the  characteristics  of  the  uncompensated  open-loop 
system  under  consideration  make  the  existence  of  such  a  compensator 
impossible. 

In  order  to  continue  the  investigation  of  the  non-linear 

problem,  the  actuator,  vane  and  aerodynamic  transfer  functions  were 

arbitrarily  changed  to  make  a  satisfactory  compromise  possible. 

The  actuator  and  vane  transfer  functions  were  set  equal  to  one,  and 

the  aerodynamic  damping  ratio   5,  was  changed  from  .0685  to  .3. 

The  block  diagram  of  the  simplified  system  is  shown  in  Fig.  17. 

The  compensator  decided  upon  was  a  lead  network-amplifier  with 

transfer  function: 

^25_ 

-i 


c(s)  =    1  +  «S 


1  +  ^  « 


K_  K  was  set  at  10  db,  again  giving  a  fairly  high  value  of 
closed-loop  peak  gain,  and  requiring  that  K~  =  2.38  db  for  unity 
closed-loop  steady-state  gain.   The  closed-loop  frequency  response 
is  plotted  in  Fig.  18.   It  has  a  satisfactory  frequency  spread. 
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SECTION  6;  THE  TRANSIENT  RESPONSE  OF  THE  NON-LINEAjk.  CLOSSD- 
LOOP  SYSTEM 

In  the  simplified  system,  that  is  for 

W(s)  =  V(s)  =  1 

Eqs.  5-1  and  5-7  become,  respectively, 

E(t)  =  Kj  o<D(t)  -  <K(t)  (6-1) 

and 

Y[   (t)  =  K2  C(D)  E(t)  (6-2) 

Equations  6-1  and  6-2  are  applicable  to  the  non-linear  system  as 
well  as  to  the  linear  one. 

Using  the  same  notation  as  in  Section  2,  the  differential 
equation  for  the  aerodynamics  is: 

T(D)  <X  =  U(D)  H  -(3o<3  (6-3) 


Where 


(6-4) 


T(D)   =  r  +  2  5,<D+   tti 

1    1       ^1 

U(D)   =  Kx  U^  (1  +  ZD)  (6-5) 

Substituting  (6-l)  in  (6-2)  gives 

»|  (t)   =  K2  C(D)   [  K5  ©<D(t)  -  o<(t)]     (6-6) 

Substituting  (6-6)  in  (6-3), 

T(D)  o<  =  K2  C(D)  U(D)  [l^  <*d  -  o<^  -  ft^  (6-7) 
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This  may  be  rearranged  in  the  form: 
[t(D)  +  K2  C(D)  U(D)]  o<  =  K2  K^  C(D)  U(d)  °<D  -  /9  o<3       (6-8) 

Now  an  approximation  to  the  solution  of  Eq.  6-8  is  the  solution 
^(t)  of: 

[T(D)  +  K2  C(D)  U(D)]^1  »  Kg  Kj  C(D)  U(d)<*d  (6-9) 

In  transfer  function  form,  Eq.  6-9  is 

^t  K~  K_  C(s)  U(s) 


^D        T(s)  +  K2  C(s)  U(s) 


L 


m. 


>  S  c(s)  : 


1  +  K0  C(s) 


U(s. 


(6-10) 


^2  ^aj     T(sJ 
Comparing  Eqs.  6-4  and  6-5  with  Eq.  5-6  we  see  that 


Hfi  ■    %*<•>  (6-n) 


Therefore: 


<=<!  K,   Kp  K,  C(s)  P(s) 

=H  S )  =  2     2    ? (6-12) 

^D  1  +  Kx  K2  C(s)  P(s) 

The  right  hand  side  of  (6-12)  is  simply  the  closed-loop  transfer 
function  of  the  linear  system  (cf.  Eq.  5-9),  so  that  °<-,(t),  the 
first  approximation  to  the  transient  response  of  the  non-linear 
system,  is  the  response  of  the  linear  system  -  not  an  unexpected 
result. 
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Now  by  analogy  with  the  argument  presented  in  Section  2,  we  take 
the  second  approximation  to  be : 


g(   (t)  =  R  ©^(t)  (6-13) 


where 


E  =   (0<^SH  (6-14) 

(°VgW  is  the  steady-state  output  of  the  non- 
linear system, 

(°()    is  the  steady-state  output  of  the  linear 

system  (equal  to  the  magnitude  of  the  input 
in  the  present  example). 

The  third  approximation  to  ©<(t)  is  °C(t),  the  solution 


of: 


[T(D)  +  K2  C(D)  U(D)]o<5  =  K2  K^  C(D)  U(D)  °<d  -  ft  <*£ 

Or,   if    ^..At)     is  the  solution  of: 
3* 

[T(D)  +  K2  C(D)  U(D)]c<^    =     -fi°<2  (6-16) 

then  by  superposition  of  the  solutions  of  (6-9)  and  (6-16), 

o<3   =  <X1     +    o^3#  (6-17) 

The  solution  of  Eq.   6-16  follows  the  same  lines  as  that  of 
Eq.   2-16:     we  have    o<,  (t)     in  the  form  of  a  Fourier  series: 

<*l(t)    =  |«rfD|  ^    +      1   ^  t*  sin  (n%t  +  Bn)] 
where     \°^ri  |       is  the  magnitude  of  the  step  input,   and  the  A's 
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and  B's  represent  the  frequency  response  of  the  closed-loop  system. 
As  before  we  let 


A  2    V    A 

o<r  (t)  =   r  +  fr  z  tt  sin  ^V  +  Bn 


) 


lo 


so  that 


(t)     =       (o<      |       o<       V 

1  1M  Ay 

From  Eq.   6-13: 

<*2(t)     =     R^lt)     =     Rlo<D|     c*     (t) 

and 

°^23    =  r51°<D|3  <<i  (6-18) 

We  find  the  Fourier  series  for  °^1  by  the  method  of  Section  2, 

and  then,  by  Eq.  6-18,  the  series  for  ©<  J .       By  Eq.  6-16,  <=<  _„ 

2  3* 

is  found  by  applying  the  frequency  response 

H(icO)  =  ^ (6-19) 

T(i«A>)  +  K  C(i*o)  U(i«j) 

to  the  series  for  -°*?  •   H(iuj)  may  be  rearranged,  using 
Eq.  6-11,  to  simplify  its  calculation.   It  becomes: 

/3  r  Kx  K2  K^  C(ito)  P(i«o)  1 

"(iw)  U(i«)  L  l  +  I^Kg  C(iu>)  P(ia>)  ]    (6"20) 


£ 


K^CO 


The  bracketed  term  is  the  linear  system  closed-loop  frequency 
response. 

The  next  approximation  to  <K(t),  if  it  is  desired  to  go 
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further,    is     °^,(t),    the  solution  of: 

^T(D)  +  K2  C(D)  U(D)]o<4  =  K2  K^  C(D)  U(D)c<D  -  ft  «<  '  (6-21 ) 

Adding  and  subtracting  p  o*' ^  on  the  right  hand  side  of  Eq.  6-21: 

£u(D)  +  K2  C(D)  U(D)]  o<4  =  K2  i^  C(DJ  U(L)  o<^ 

-poC23+/3(«23-o<33)        (6-22) 
Let   °^,^(t)  be  the  solution  of: 

[T(D)  +  K2  C(D)  U(D)]  o<^  =  /3  (  «^3  _  0^3)  (6_23) 

Then  by  superposition  of  the  solutions  of  (6-15)  and  (6-23), 

To  solve  Eq.  6-23,  the  Fourier  series  for  (  °(        -  c<   )  is 
found  by  numerical  harmonic  analysis  in  exactly  the  same  way  as 
outlined  in  Section  2  and  demonstrated  in  Section  3.   The 
frequency  response  H(icu)  is  applied  to  it  to  give  °^A^>  and 
c><   is  found  by  adding  the  evaluations  of  c<   and  <=><.,*• 
Further  iterations  follow  the  same  lines. 

Having  calculated  °^-,(t)  and  °<  #(t)  for  a  step  function 
of  any  magnitude,  it  is  easy  to  find  their  values  for  any  other 
input  magnitude.   Let  the  subscript  "a"  denote  the  case  which 
has  been  calculated  and  "b"  the  case  which  is  to  be  calculated. 
The  relationship  between  the  linear  solutions  is,  of  course, 

o<  (t)  =   ^Dbl    o<1  (t)  (6-25) 
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From  Eq.   6-13, 

o<2   (t)       =     Ra    o<1   (t)  (6-26) 


a 

and 


°<2   (t)       =     Rb    <^±   (t)  (6-27) 

b  b 

Substituting  (6-25)  in  (6-27): 

«2   (t)       =     \  '!^b[      ^   (t)  (6-28) 

l°\  I  a 

Eliminating  o<  (t)  between  (6-26)  and  (6-28): 
a 

*<2  (t)   =   7)l   Db'    o<2  (t)  (6-29) 

Thus 


R  ' 
a 

a 


o<2-     =   u>  "bl\      ^23  (6-30) 


Hd, 


R 
a 

a 


Due  to  the  linearity  of  Eq.  6-16, 


,h  Ki\3 


°<3.     =     '    b'\    <*,»  (6-51) 

Having  found  o<   (t)  and  c<_#  (t),   c*   (t)  follows  from 
••■b  ->  b        ->b 

Eq.  6-17.   No  new  Fourier  evaluations  have  been  performed. 

The  only  simplification  possible  in  the  next  iteration  is 

3  3 

that  the  evaluation  of  °^2v   follows  from  Eq.  6-13.   ©<  * 

must  be  found  by  cubing  o<    point  by  point,  and  a  new  harmonic 

•?b 
analysis  must  be  performed. 
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We  now  apply  this  procedure  to  the  simplified  control 
system,  first  for  a  step  input  °<   of  magnitude  .3  radian. 
The  fundamental  frequency  u>p  is  taken  as  .38  cO^,  one  fifth 
of  the  closed-loop  frequency  response  peak  in  Fig.  18.   The 
tabulations  of  Fourier  series  and  frequency  responses  for  this 
example  will  be  found  in  Table  XIII,  and  the  evaluations  in 
Table  IXV.   The  Fourier  series  for  a  square  wave,  magnitude  .3, 
is  first  tabulated,  then  the  closed-loop  frequency  response. 
The  Fourier  series  for  the  output  of  the  linear  system,  °^-|(t), 
is  found  by  applying  the  closed-loop  frequency  response  to  the 
square  wave.   We  find  the  steady-state  output  of  the  non-linear 
system  from  Eq.  5-H>  which,  upon  substitution  of  the  constants, 
becomes: 

(*),M+.56  («)'      =     |^,|  (6-32) 


For 


find 


(°<>SL  -   Kl=  •' 


(*)SN    =       .287 


Thus,   from  Eq.   6-14, 


R    =      -^      =     .957 
.3 


°<p(t),  found  by  Eq.  6-13,  is  evaluated. 

The  Fourier  series  for  °C   ,  found  by  the  method  developed 

1  u 
in  Section  2,  is  tabulated. 

r3  I**])!3  -  (-957)5  (.3)3  =  .02365 

By  Eq.  6-18,  -  °<3  is  found  from  °<^    and  tabulated.   The 
2  1  U 

frequency  response  H(ito)  is  tabulated,  and  the  Fourier  series 
for   °C*  is  obtained  by  Eq.  6-16,  tabulated,  and  evaluated. 
By  Eq.  6-17,  the  evaluation  for  °C(t)  is  found  by  adding  those 
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for  °<,   and   «*_. 

1  3* 

°<2(t)  and  ©C(t)  are  plotted  in  Fig.  19.   We  note  that 
c*  (t)  is  not  very  different  from  °C(t)  and  is  therefore 
probably  a  close  approximation  to  <=<(t).   The  effect  of  the 
non-linearity  is  small,  causing-  a  very  slight  increase  in  frequency. 

To  illustrate  a  case  where  the  effect  of  the  non-linearity 
is  large,  we  take   (^rJ   =  1  radian.   Although  this  is  an 
unrealistic  value  of  angle-of-attack,  the  result  will  be  equivalent 
to  that  of  a  missile  at  low  angle-of-attack  but  with  a  larger 
non-linearity.   It  is  simply  convenient  to  calculate  it  in  this 
manner  since  the  equations  have  been  formulated  on  this  basis. 
The  evaluations  are  given  in  Table  IXV.   From  Eq.  6-32  find 


(-OgK   =  .757 


Using  the  subscript  "b"  to  denote  this  case,  we  have,  from 
Eq.  6-14, 

\    =   -^P   =  .757 
The  subscript  "a"  denotes  the  .3  input  case,  already  evaluated. 

b    ^       (.757)  (1)       2fi- 
Ra  Ka|  =    trntS)       ~       2'63 

The  evaluation  of  oC0   (t),  found  according  to  Eq.  6-29,  is 

^b 
tabulated . 

S    l~I>al/ 

©C*  (t)  is  found  by  Eq.  6-31,  and  ©C   by  Eq.  6-17.    °<2  and 

°^3b  are  plotted  in  Fig.  20. 
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One  additional  graph  is  presented  in  Fig.  21.   Here  the  third 
approximations  corresponding  to  both  the  .3  and  the  1.0  input  have 
been  plotted,  as  has  that  for  the  .65  case,  the  calculations  for 
which  are  not  shown.   These  solutions  have  been  "normalised'1  by- 
dividing  each  ordinate  by  the  corresponding  steady-state  value 
(o<)      °^p(t)/(°<)„AI  ,  which,  of  course,  is  the  same  for  all 
inputs,  has  also  been  plotted.   This  is  a  useful  method  for 
visualising  the  effect  of  increasing  input  magnitude  on  °C(t). 
It  would  also  be  useful  for  plotting  other  approximations,  or  the 
exact  solutions. 

Improved  approximations  will  not  be  demonstrated  here,  and 
the  determination  of  the  number  of  approximations  necessary  for 
acceptable  convergence  remains  for  the  worker  who  thinks  this  method 
might  be  useful. 

Finally,  it  is  interesting  to  note  that  a  linear  closed-loop 
system  with  two  inputs  can  be  described  which  is  equivalent  to  the 
non-linear  system.   The  block  diagram  of  this  system  is  shown  in 
Fig.  22.   Comparison  with  Eq.  6-8  will  prove  its  equivalence. 

Fig.  22  is,  in  fact,  the  general  outline  of  a  possible  analog 

3 

computer  set-up  where     c<       is  found  by  putting     oC  through  a  cubing 

device. 


50 


DISCUSSION 

It  has  been  mentioned  that  a  quick  method  for  determining  the 
effect  of  non-linear  stiffness  was  desired.   Specifically,  it  was 
hoped  that  few  terms  would  be  required  in  the  first  iteration,  and 
that  one  iteration  would  be  sufficient  for  acceptable  accuracy,  at 
least  in  the  first  overshoot  region.   In  fact,  the  number  of 
significant  terms  was  not  small,  and  a  single  iteration  was  not 
sufficient  no  matter  how  many  terms  were  retained,  except  for  small 
inputs. 

The  frequency  spread  requirement,  discussed  in  Appendix  II, 
places  a  limitation  on  the  range  of  systems  to  which  the  Wass, 
Hayman  method  may  be  applied.   And  the  Wass,  Hayman  method  must  be 
applicable  for  the  present  method  to  work.   If  the  frequency  spread 
is  of  the  order  of  20,  the  number  of  terms  in  the  Fourier  series 
could  be  increased  to  handle  it,  but  if  the  original  control  system 
designed  here  is  at  all  typical,  the  frequency  spread  is  likely  to 
be  at  least  50,  requiring  at  least  25  terms  in  the  Fourier  series 
for  the  linear  response,  and  more  for  the  non-linear  response. 
The  need  would  have  to  be  great  to  justify  the  effort  of  evaluating 
such  series. 

In  Section  5  an  unrealistic  value  of  aerodynamic  damping  was 
arbitrarily  chosen  to  produce  a  control  system  to  wnich  the  present 
method  could  be  applied.   It  is  possible  that  a  more  justifiable 
change  in  the  control  system  could  have  been  made,  but  this  was  not 
done  due  to  a  necessity  to  find  an  acceptable  system  quickly.   It 
remains  to  be  seen  whether  the  present  method  can  be  used  in  any 
practical  cases. 

Although  only  "hard"  stiffness  has  been  considered,  i.e. 
positive  /3  ,  the  method  is  equally  applicable  to  soft  stiffness, 
or  negative  ft   .   In  the  latter  case  the  presence  of  static 
instability  at  high  angle-of -attack  must  be  considered. 


51 


CONCLUSIONS 


1.  Use  of  the  iterative  technique  developed  is 
capable  of  giving  the  transient  response,  to  an  acceptable  degree 
of  accuracy,  if  the  frequency  response  of  the  system  meets  the 
frequency  spread  requirement. 

2.  In  an  example,  acceptable  accuracy  was  obtained 
after  two  iterations.   The  number  of  iterations  necessary  when 
the  non-linear  effect  is  large  was  not  determined. 

3.  The  labor  required  to  evaluate  the  Fourier  series 
and  to  perform  the  numerical  harmonic  analyses  necessary  in  the 
second  and  subsequent  iterations  is  considerable,  thus  limiting 
the  practical  value  of  the  method. 
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APl;EHDIX  I 

CALCULATION  OF  THE  AERODYNAMIC  CONSTaKTS 

Brown-Edwards,  in  Ref .  1,  found,  for  the  missile  under 
consideration  and  for  a  e.g.  position  8  body  diameters  aft  of  the 
body  apex,  the  following  values: 

CL    =     2.87  *  -  .222  n  (1-1) 

CM    =  -  1.14  <*  -  2.65  <*3  +  .97  n         (1-2) 
eg. 

Fixing  the  body  diameter  1  at  one  foot,  he  estimated  for  a 
typical  missile: 

^  to  ak  lb. sec 

rtn,  =   12.05  — ji — 

and 

I  =  250  lb. ft. sec.2 

2 
The  reference  area  S  IS   9  f t  . 

We  take  the  case  of  flight  at  50,000'  at  Maeh  3.   The  air 

2  /  4 
density  p      is  .000362  lb. sec  /ft   and  the  speed  of  sound  c 

is  968  ft/sec.   Thus  the  velocity  is: 

V  =  3  (968)  =  2900  ft/sec. 
The  dynamic  pressure  is 

Q  =  ^V2  =  i( •  000362 ) (2900)2  =  1520  3 


ft' 


and 


QS  =  (1320) (9)  =  13,700  lb. 


Now 

2     *>   V1    =  "  d  s    ^3    =     -  (13,700;(2.87)  =     -  39,400  lb. 

where       "3C  /3c<  is  obtained  from  Eq.    1-1. 

2 |ii=  -    QS     r^     =     -   U3,700)(-   .222)  =  3,040  lb. 

where      ()C  /<)n  is  obtained  from  Eq.    1-1 

Mn     =       |I  -  qsI  — ^'  =  (13,700)(1)(.97)  =  13,300  lb. ft. 


where  3C     /g,.   is  obtained  from  Eq.  1-2.   The  moment  due 
e.g.   * 

to  ©<  is,  from  Eq.  1-2, 

M"(<*)   =  Q^2     (-1.14  <*  -  2.65  <*?) 

=   (13,700)(1)(-  1.14<X-  2.65<*5) 

(-15,600  oC  -  36,300©C5)  lb. ft. 
Comparing  this  with  Eq.   1-14,   we  see  that 
M«       =     -15,600  lb. ft. 


and 


M  =    -36,300  lb.ft. 


Using  Eqs.   1-17  to  1-21, 


"^i2       =     -      ^?    =       15,600      =     67.8  sec"  , 
I  230 


"S 

= 

8.23  sec"1 

*1 

= 

£« 

39,400 

1 

2(12.05)(2900)(8.23) 

= 

.0685 

h 

= 

M 

1 

13,300            =       .853 
(230)(67.8) 

z 

= 

Zn 

3,040 

<y*  V  K,    ^ 

1            1 

(12.05)(2900)(.853)(67.8) 

= 

.001505  sec. 

ft 

= 

M*3 
I 

^2SSL    =         15Q  sec-2 
230 

APPENDIX  II 

THE  WASS,  HAYMAN  METHOD  OF  DERIVING  THE  TRANSIENT  RESPONSE  OF 

A  LINEAR  SYSTEM  FROM  THE  FREQUENCY  RESPONSE 

Reference  2  presents  a  method  of  deriving  the  approximate 
response  of  a  linear  system  to  a  step  input  from  the  frequency- 
response.   The  method  is  to  find  the  transient  response  to  a 
square  wave  input  instead  of  to  a  step  input.   The  justification 
for  such  a  manouver  is  best  explained  with  reference  to  Fig.  23  . 
Figure  23a  is  the  graph  of  a  step  input,  and  Fig.  23b  is  the 
corresponding  output  of  a  typical  system.   Figure  23c  is  a  square 
wave  input  of  period  T,  and  Fig.  23d  is  the  corresponding  output. 
It  is  evident  that  if  the  half -period  T/2  is  greater  than  the 
time  required  for  the  transient  to  completely  damp  out,  then  any 
of  the  half  cycles  of  Fig.  23d  is  the  same  as  the  transient  part 
of  Fig.  23b.   Of  course  in  linear  systems  the  transient  never 
completely  damps  out.   Therefore  Fig.  23d  is  really  an 
approximation  to  Fig.  23b,  its  accuracy  depending  on  the  extent 
to  which  the  transient  subsides  in  the  time  T/2. 

The  reason  for  using  a  square  wave  is  that  it  can  be  expressed 

as  a  Fourier  series.   The  Fourier  series  for  a  square  wave  F.(t) 

with  equal  mark-to-space  ratio,  minimum  value  zero,  maximum  value 

unity,  and  fundamental  frequency  «J  ,  can  be  written: 

a 


Fi(t)  =  i+;    2 


1  sin   n  i*Lt 

n         P 


n  =  1,  3,  5,  ...,  «*» 

(Wass  and  Hayman  use  the  symbol  n  differently  -  this  will  be 
explained  shortly).   If  F±(t)  is  the  input  to  a  linear  system 
the  magnitude  and  phase  of  each  component  will  be  changed 
according  to  the  frequency  response  of  the  system,  and  the 


output  will  be 


=    —  +   -   Z  —  am  (n<Opt  +  Bn) 
n  =  1,  3,  5,  ...  <*> 


where   A   is  the  gain  of  the  system  at  frequency   n  tu  , 

n  Jr 

B   is  the  phase  shift  of  the  system  at  frequency  ntO  . 

In  order  to  use  this  method  numerically,  the  summation  must 
be  limited  to  a  reasonably  low  number  of  terms.   Wass  and  Hayman 
state  that  it  has  been  found  by  experience  that  a  satisfactory 
compromise  between  labor  and  accuracy  is  achieved  by  considering 
up  to  the  11th,  or  sometimes  13th,  harmonic  of  the  fundamental 
square  wave  frequency.   Having  decided  on  the  number  of  terms 
to  be  summed,  a  value  for  the  fundamental  frequency  to   must  be 
chosen.   Wass  and  Hayman  suggest  that  for  systems  whose  gain 

curves  exhibit  a  definite  peak,  tu  should  be  one  fifth  of  the 

J? 

frequency  of  the  peak.   Suggestions  are  also  made  with  respect 
to  other  types  of  systems. 

mass  and  Hayman  do  not  point  out  that  taking  the  summations 
up  to  even  the  13th  harmonic  will  fail  to  give  satisfactory 
accuracy  in  certain  types  of  systems.   They  do  infer  that  the 
transient  response  is  mainly  determined  by  the  characteristics  of 
the  frequency  response  at  frequencies  lower  than  that  at  which 
the  gain  is  15  db  below  the  zero  frequency  gain.   They  do  not 
point  out  that  the  gain  and  phase  shift  at  the  fundamental 
frequency  must  not  differ  greatly  from  their  zero  frequency  values. 
For  conciseness  we  may  use  the  term  "frequency  spread"  to  denote 
the  ratio  of  the  -15  db  frequency  to  the  frequency  at  which  the 
gain  or  phase  first  differs  appreciably  from  its  zero  frequency 
value.   Thus  the  Wass,  Hayman  method,  summing  up  to  the  13th 


harmonic,  is  not  applicable  to  systems  whose  frequency  spread  is 
greater  than  13. 

An  important  part  of  the  Wass,  Hayman  report  is  a  table,  the 
"if"  Table,  which  may  be  used  to  simplify  the  evaluation  of  Fourier 
series.   This  is  a  table  giving  the  values  of  mo  t'  in  degrees, 
for  39  values  of  non-dimensional  time  t',  and  for  odd  values  of 
n  from  1  to  13.   Even  values  of  n  are  not  included  because  the 
Fourier  series  for  the  responses  of  linear  systems  to  square  wave 
inputs  contain  only  odd  harmonics.   Since  the  evaluation  of  even 
harmonics  is  required  in  this  thesis,  an  extension  to  the  if  Table 
is  required,  and  this  is  given,  for  even  values  of  n  from  2  to  8, 
in  Table  XV.   Wass  and  Hayman' s  use  of  the  symbol  n  precludes 
the  possibility  of  even  harmonics,  and  this  is  the  reason  it  is 
used  differently  here. 

Wass  and  Hayman  present  a  design  for  a  circular  computer  which 
reduces  considerably  the  time  required  to  evaluate  a  Fourier 
series.   The  computer  suggested  consists  of  a  circular  base  with 
various  lines  and  scales  marked  on  it,  two  arms  which  pivot  at  the 
center  of  the  circular  base,  and  a  cursor  which  slides  along  one 
of  the  arms.   The  author  constructed  and  used  such  a  computer,  and 
eventually  devised  a  different  design  which  was  superior  in  several 
respects.   Fig.  24  is  a  photograph  of  this  computer.   It  was 
originally  a  "Wind  and  Navigation  Calculator  Mk.  I",  and  all  the 
scales  around  the  periphery  are  superfluous  except  the  360  degree 
one.   It  consists  of  an  opaque  base  and  a  transparent  top, 
pivoted  at  the  center.   The  horizontal  lines,  the  numbers  running 
down  the  center,  and  the  arrowhead  are  etched  on  the  base.   The 
360  scale  is  etched  on  the  top.   The  top  is  roughened  so  that  it 
will  take  pencil  marks.   To  evaluate  the  term 

A  =  3  sin  (mo.-.f  +  ©  ) 

if 

the  procedure  is  as  follows: 


1.  Rotate  the  top  until  the  arrowhead  points  to 
(90°  -  0  )  on  the  360°  scale. 

2.  Make  a  pencil  mark,  on  the  top,  at  the  point 
where  the  value  of  B  appears  on  one  of  the  vertical  scales 
above  the  center. 

3.  Rotate  the  top  until  the  arrowhead  points  to 
nii/pf  on  the  360  scale. 

4.  Follow  the  new  position  of  the  pencil  mark  along 
the  horizontal  lines  to  the  same  scale  used  to  place  the  mark, 
and  read  A.   It  is  positive  if  above  center,  negative  if  below. 

5.  Repeat  steps  3  and  4  for  other  values  of  nu>  t1. 

The  new  computer  is  superior  to  the  Wass,  Hayman  version  in 
these  respects: 

a.  The  two  arms  of  the  Wass,  Hayman  computer  must  be  turned 
through  several  revolutions  without  moving  relative  to  each  other, 
yet  the  angle  between  them  must  be  adjustable.   A  design 
incorporating  such  a  feature  may  be  difficult  to  achieve,  particularly 
for  the  home  builder.   The  new  computer  obviates  this  difficulty. 

b.  The  angle*  between  the  arms  and  the  position  of  the 
cursor  must  be  reset  for  each  new  term  of  the  Fourier  series. 
With  the  new  computer  a  mark  for  each  term  of  the  series  can  be 
made  at  the  outset,  each  one  being  labeled  appropriately.   This 
may  not  seem  important,  but  in  practice  it  is.   Mistakes  can  be 
corrected  easily  because  the  setting  is  not  "erased"  for  the  next 
term  as  with  the  Wass,  Hayman  version,  and  evaluations  can  be  made 
for  additional  values  of  time,  without  resetting,  after  several 
points  have  been  plotted. 


APPENDIX     III 
GRiiPHICAL  ADDITION  Qg  LIKE  riABHQNICo 

If,    in  the  expression 

e  sin  (u>t  +  ©)  +  f  sin  («*>t  +   <¥>)  +  g  sin  (u^t  +  if  ) 

e,  f,  g,  G,     <p ,   and  vp  are  known  constants,  then  an  equivalent 
expression  is 

h  sin  ( u)  t  +  ri  ) 

where  h  and  »j  are  determined  according  to  the  construction 
below. 


This  method  is  valid  for  any  number  of  terms  in  the  original 
expression. 


APPENDIX  IV 

TRANSFORATION  OF  THE  CUBIC  EQUATION 
Given  the  equation: 


d2%  d* 

dt2        1    dt 


+  ao^+  bX;  =   q       (i-i) 


Make  the  substitutions: 


t 

1 
fa 

1      0 

-  z 

and 

X 

■V: 

-   Is 

Then 

Eq. 

1-1  becomes: 

<LZ2 

+       al 

0 

<LZ 

We  shall 

call  the 

quantity 

\ 

Q 

rs 

a  54 

0 

+   3        Qfb 
2   ^         "T7/2 


o 


a 


the  non-dimensional  input.   The  coefficient   ==-   is  equal  to 

o 

twice  the  damping  ratio  5  • 


APPENDIX  V 

DESIGN  OF  TiJE  INCIDENCE  VANE 

The  vane  is  conical,  and  it  is  necessary  that  the  e.g.  be  as 
far  forward  as  possible.   The  configuration  to  be  used  is  that 
shown  in  Fig.  25.   The  nose  will  be  solid  steel,  and  the  base  will 
be  a  thin  shell  of  light-weight  plastic.   The  optimum  ratio  of 
nose  length  to  tail  length,  for  most  forward  e.g.  position,  is  that 
which  places  the  e.g.  at  the  junction  of  the  nose  and  tail.   The 
equation  for  this  condition  was  found  to  be: 


i   Pi 


t     1   '  1        sin 


1  4  /°2    (2n4  -  3„3  +  n) 

where        t  is  the  base  thickness, 

L  is  the  cone  length, 

f*-,  is  the  density  of  the  nose  material, 

/>  is  the  density  of  the  base  material, 

©  is  the  cone  semi-angle, 

—   is  the  distance  of  the  e.g.  from  the  apex. 

The  density  of  steel  is  about  .3  lb/cu  in.,  and  that  of  a  high- 
temperature  silicone  plastic  is  about  .05  lb/cu  in.  (p.  66  of 
Ref.  4).   Using  these  values,  the  compromise  chosen  was: 

0  =  22.5° 

n  =  3 

L  =  3" 

t  =  . 0614" 

It  is  questionable  whether  the  base  is  thick  enough  for  adequate 
strength,  but  it  will  be  assumed  that  such  a  design  is  practicable, 


APPENDIX  VI 


CALCULATION  OF  THE  VANE  CONSTANTS 

Ref .  5  gives  experimental  pitching  moment  data  for  cones  in 
supersonic  flow  in  the  form  of  graphs  of  C.,  versus  °<  ,  where 
the  pitching  moment  about  the  nose  is 

^OSE  =   ^^Tllc°ll 


and 

r,   =  base  radius, 

1   =  cone  length. 


The  curves  of  CL  versus  <*  were  virtually  linear  for  the  conditions 
applicable  to  the  present  case,  and  the  value  of  "9CT./3°<   about 
tne  nose  was  .020.   The  center  of  pressure  remained  very  close  to  a 
position  2/3  1   back  from  the  nose,  so  that  the  moment  curve  slope 
for  a  cone  hinged  at  l/j  1  would  be  half  of  .020,  or  .010. 

From  Appendix  I: 


Q 

=  1520  lb/ft2 

From  Appendix  v 

1 
c 

=  .25  ft. 

< 

=  .0107  ft2. 

Thus : 


/2M 

1     la^/ 


I-  7)  n-Q/r2 


a -j 


i'i 


=  -  TT (1520) (.0107) (.25) (.01) 

=  -  .128  ft. lb. 
Iy    =     1.70  x  10"5  ft. lb. sec? 

(Found  by  integration,  calculation  not  shown). 


uJn2         =       -    fl       =  -128 =       7,550  sec'' 

2  h  1.70    x     10"5 


U>  =      86.9  sec"1 

52        =      .5 

(Chosen  arbitrarily  on  the  assumption  that  the  vane  damping  device 
can  be  so  designed. ) 


C2     =     2Ivu>n^2     =     2(1.70  x  10"5)(86.9)(. 5) 


=     .00148  ft. lb. sec. 


Native  method  for  finding  the  tran 
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